#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>

#include "Vector/vector.03.h"
#include "Color/color.01.h"
#include "Timer/timer.01.h"
#include <vector>
#include <deque> 

#define SQR(x) ((x)*(x))

using namespace std;

typedef unsigned int uint;
typedef long unsigned int luint;
typedef double real;



const double kB = 1.3806488e-23;

const int gsizex = 50;
const int gsizey = 50;
const int gsizez = 50;

double gstepx = 3;
double &gstepy = gstepx;
double &gstepz = gstepx;



double drij2max = SQR(3.0); 
double dt = 1e-3;
double period = 15;
uint nsteps = period/dt;

double lU = 1e-9; // m [1nm]
double eU = 1.60217646e-19; // J [1ev]
double mU = 1.660538921e-27; // kg [amu]
double tU = lU*sqrt(mU/eU); // sec

luint memmax = 500; //[MB]	

bool ssF = false;
double sstime = 0.01;

bool coolingF = false;

double all_memsize;

int ii(1), e00(1), e11(1), e01(1);


struct P {
	Vector<real> R;
	Vector<real> V;
	Vector<real> F;
	int type;
};
vector<P> p;
vector<P*> g[2*gsizex][2*gsizey][2*gsizez];

struct PW {
	double e;
	double s; 
	double rc;
};
PW pw[10][10];

struct PT {		
	double m;
	Color color;	
};
vector< PT > pt;		

uint ntot;
				
Timer tr, trf;
		
double rnd() {
	return (double)rand()/RAND_MAX;
};


template <int ROW, int COL>
void packToHCP(vector<P> &p, const PW (&pw)[ROW][COL], int type,  Vector<double> R, Vector<int> S, Vector<double> V, double t) {
			
	double dr = pw[type][type].s*0.977;
			
	for (int i=0; i<S.x; i++)
		for (int j=0; j<S.y; j++) 
			for (int k=0; k<S.z; k++){
					
				P pi;	
					
				pi.R.x = dr*i + (j%2==0 ? 0: 0.5*dr) + (k%2==0 ? 0: 0.5*dr);
				pi.R.y = 0.5*sqrt(3)*dr*j + (k%2==0 ? 0: sqrt(3.0)/6.0*dr);
				pi.R.z = sqrt(6.0)/3.0*dr*k;
					
				pi.V = Vector<double>(t*(1-2*rnd()), t*(1-2*rnd()), t*(1-2*rnd()));

				pi.R  += R;
						
				pi.V  += V;
				pi.type = type;	
				p.push_back(pi);	
					
			}
}


template <int ROW, int COL>
void packToHCP(vector<P> &p, const PW (&pw)[ROW][COL], int type1, int type2, Vector<double> R, Vector<int> S, Vector<double> V, double t) {
			
	double dr = pw[type1][type1].s*0.977;
			
	for (int i=0; i<S.x; i++)
		for (int j=0; j<S.y; j++) 
			for (int k=0; k<S.z; k++){
					
				P pi;	
					
				pi.R.x = dr*i + (j%2==0 ? 0: 0.5*dr) + (k%2==0 ? 0: 0.5*dr);
				pi.R.y = 0.5*sqrt(3)*dr*j + (k%2==0 ? 0: sqrt(3.0)/6.0*dr);
				pi.R.z = sqrt(6.0)/3.0*dr*k;
					
				pi.V = Vector<double>(t*(1-2*rnd()), t*(1-2*rnd()), t*(1-2*rnd()));
				
				pi.R.x-= 0.5*S.x*dr;
				pi.R.y-= 0.5*S.y*0.5*sqrt(3)*dr;
				pi.R.z-= 0.5*S.z*sqrt(6.0)/3.0*dr;
				
				pi.R  += R;
						
				pi.V  += V;
				
				pi.type =  (i==0 || j==0 || k==0 || i==S.x-1 || j==S.y-1 || k==S.z-1) ? type2 : type1;
				
				p.push_back(pi);	
					
			}
}


void calc_force(vector<P*> &pi) {
	
	
	for (int i=0;i<pi.size()-1;i++)
		for (int j=i+1;j<pi.size();j++) {
						
			real dx = pi[i]->R.x - pi[j]->R.x; // 1
			real dy = pi[i]->R.y - pi[j]->R.y; // 1
			real dz = pi[i]->R.z - pi[j]->R.z; // 1
								
			real r2 = dx*dx + dy*dy + dz*dz; // 5

			if ( r2 > pw[pi[i]->type][pi[j]->type].rc ) continue; 
							
			real r4 	= r2*r2; // 1
			real r6i 	= 1.0/(r2*r4); // 2
			real r12i = r6i*r6i;  // 1
			
			real f = pw[pi[i]->type][pi[j]->type].e*r6i*r4*(r12i-r6i); // 4
							
			real dxf = dx*f;
			real dyf = dy*f;
			real dzf = dz*f;
							
			pi[i]->F.x += dxf; // 2
			pi[i]->F.y += dyf; // 2
			pi[i]->F.z += dzf; // 2
						
			pi[j]->F.x -= dxf; // 2
			pi[j]->F.y -= dyf; // 2
			pi[j]->F.z -= dzf; // 2
							
	}
}


void calc_force(vector<P*> &pi,  vector<P*> &pj) {
					
	for (int i=0;i<pi.size();i++)
		for (int j=0;j<pj.size();j++) {
						
						
			real dx = pi[i]->R.x - pj[j]->R.x; // 1
			real dy = pi[i]->R.y - pj[j]->R.y; // 1
			real dz = pi[i]->R.z - pj[j]->R.z; // 1
								
			real r2 = dx*dx + dy*dy + dz*dz; // 5

			if ( r2 > pw[pi[i]->type][pj[j]->type].rc ) continue; 
							
			real r4 	= r2*r2; // 1
			real r6i  = 1.0/(r2*r4); // 2
			real r12i = r6i*r6i;  // 1
			
			real f = pw[pi[i]->type][pj[j]->type].e*r6i*r4*(r12i-r6i); // 4
								
			pi[i]->F.x += dx*f; // 2
			pi[i]->F.y += dy*f; // 2
			pi[i]->F.z += dz*f; // 2
							
	}
					
}


void run(){
					
	for(int i=0; i<p.size(); i++) p[i].F*=0;

	uint step(0);
	
	double glenx = gsizex*gstepx;
	double gleny = gsizey*gstepy;
	double glenz = gsizez*gstepz;
	
	double g2sizex = 2*gsizex;
	double g2sizey = 2*gsizey;
	double g2sizez = 2*gsizez;
	
	do {

		tr.set();
						
		for(int i=0; i<p.size(); i++) {
			p[i].R+= (p[i].V + .5*p[i].F*dt)*dt;
			p[i].V+= .5*p[i].F*dt;
		}
						
		// boundaries
		for(int i=0; i<p.size(); i++) {
			
			if (p[i].R.x < -glenx) { p[i].R.x=-glenx+.0000001;  p[i].V.x*=-1; } else 
			if (+glenx < p[i].R.x) { p[i].R.x=+glenx-.0000001;  p[i].V.x*=-1; }
			
			if (p[i].R.y < -gleny) { p[i].R.y=-gleny+.0000001;  p[i].V.y*=-1; } else 
			if (+gleny < p[i].R.y) { p[i].R.y=+gleny-.0000001;  p[i].V.y*=-1; }
			
			if (p[i].R.z < -glenz) { p[i].R.z=-glenz+.0000001;  p[i].V.z*=-1; } else 
			if (+glenz < p[i].R.z) { p[i].R.z=+glenz-.0000001;  p[i].V.z*=-1; }
			
		}
						
		
		for (int xi=0; xi<g2sizex; xi++)
		for (int yi=0; yi<g2sizey; yi++)
		for (int zi=0; zi<g2sizez; zi++)
			g[xi][yi][zi].resize(0);
											
		for(int i=0; i<p.size(); i++)
			g[int(p[i].R.x/gstepx)+gsizex][int(p[i].R.y/gstepy)+gsizey][int(p[i].R.z/gstepz)+gsizez].push_back(&p[i]);
		
						
		trf.set();
						
		for (int xi=0; xi<g2sizex; xi++)
		for (int yi=0; yi<g2sizey; yi++)
		for (int zi=0; zi<g2sizez; zi++) {
							
			vector<P*> &pi = g[xi][yi][zi];
			
			if (pi.size()==0) continue;
								
			for (int i=0; i<pi.size(); i++) pi[i]->F*=0;
							
			calc_force(pi);

							
			for (int xj= (xi==0 ? xi:xi-1); xj<=(xi<g2sizex-1 ? xi+1: g2sizex-1); xj++)
			for (int yj= (yi==0 ? yi:yi-1); yj<=(yi<g2sizey-1 ? yi+1: g2sizey-1); yj++)
			for (int zj= (zi==0 ? zi:zi-1); zj<=(zi<g2sizez-1 ? zi+1: g2sizez-1); zj++) 
			if (xi!=xj || yi!=yj || zi!=zj )
				calc_force(pi,  g[xj][yj][zj]);
							
		}
		trf.get();
						
		for(int i=0; i<p.size(); i++)
			p[i].V+= .5*p[i].F*dt;
						
		tr.get();		
		
		
		
		/*
		if (ssF && step%uint(sstime/dt)==0) {
			fprintf(stderr,"ss::output\n");
			FILE* fp = fopen("tmp.txt", "w");
			for (int i=0; i<p.size(); i++)
				fprintf(fp, "%.4e %.4e %.4e %i\n", p[i].R.x, p[i].R.y, p[i].R.z, p[i].type);
				//fprintf(stdout, "%.4Le %.4Le %.4Le\n", p[i].R.x, p[i].R.y, p[i].R.z);
			fclose(fp);
		}
		*/
		
		if (step%uint(sstime/dt)==0) {
			
			char foldername[256];
			sprintf(foldername, "%i.%i.%i.%i", ii, e00, e11, e01);
			mkdir(foldername, 0777);
			
			char filename[256];
			sprintf(filename, "%s/ss.%.1f", foldername, step*dt);
			
			printf("%s\n", filename);
			
			FILE *fp = fopen (filename , "wb");
		
			struct MEM { Vector<float> R; float type;} *mem;
			mem = new MEM[p.size()];
			for (int i=0; i<p.size(); i++) {
				mem[i].R = p[i].R;
				mem[i].type = p[i].type;
			}
		
			fwrite(mem, sizeof(MEM), p.size(), fp);
			fclose(fp);
			delete [] mem;
		}
		
		
		step++;
						
		//fprintf(stderr,"%9i %10.4f  %9i %10.4f %10.4f %10.2e  %10.4f  %10.4f\n",step, step*dt,  p.size(), trf.sin, tr.sin, p.size()/tr.sin, 1/tr.sin, all_memsize);
						
	} while (step<nsteps);
	
}


int main(int argc, char *argv[]) {
	
	uint boxsize=64;
	
	
	ssF=true;
	
	period = 2.0;
	dt = 0.00025;
	nsteps=period/dt;
	gstepx = 3;
	sstime = 0.01;
	
	
	for (int i=1;i<argc;i++)
		     if (!strcmp(argv[i],"-ssF")) { ssF=true;  mkdir("ss", 0777);  }
		else if (!strcmp(argv[i],"-nsteps")) nsteps = atoi(argv[++i]);
		else if (!strcmp(argv[i],"-gstepx")) gstepx = atof(argv[++i]);
		else if (!strcmp(argv[i],"-dt"))     { dt = atof(argv[++i]); nsteps=period/dt;  }
		else if (!strcmp(argv[i],"-period")) { period=atof(argv[++i]); nsteps=period/dt; }
		else if (!strcmp(argv[i],"-boxsize")) boxsize=atoi(argv[++i]);
		else if (!strcmp(argv[i],"-sstime")) sstime=atof(argv[++i]);
		else if (!strcmp(argv[i],"-e00")) e00 = atoi(argv[++i]);
		else if (!strcmp(argv[i],"-e11")) e11 = atoi(argv[++i]);
		else if (!strcmp(argv[i],"-e01")) e01 = atoi(argv[++i]);
		else if (!strcmp(argv[i],"-ii")) ii = atoi(argv[++i]);
	
	pt.resize(2);

	
	pt[0].m = 1;
	pt[1].m = 1;
	
	pw[0][0].s  = 1;
	pw[0][0].e  = e00; // 10000
	pw[0][0].rc = 3;
	
	pw[1][1].s  = 1;
	pw[1][1].e  = e11; // 200
	pw[1][1].rc = 3;
	
	pw[0][1].s  = pw[1][0].s = 1;
	pw[0][1].e  = pw[1][0].e = e01; // 1000
	pw[0][1].rc = pw[1][0].rc = 3;	
	
	
	
	packToHCP(p, pw, 0, 1, Vector<double>(-0.5*boxsize-3,0,0), Vector<int>(boxsize,boxsize,boxsize),  Vector<double>(+15,0,0), 0);
	packToHCP(p, pw, 0, 1, Vector<double>(+0.5*boxsize+3,0,0), Vector<int>(boxsize,boxsize,boxsize),  Vector<double>(-15,0,0), 0);
	//packToHCP(p, pw, 0, 1, Vector<double>(+boxsize+2,0,0), Vector<int>(boxsize,boxsize,boxsize),  Vector<double>(-30,0,0), 0);
	//packToHCP(p, pw, 1, Vector<double>(4+boxsize+2, 4,4), Vector<int>(2,16,16),  Vector<double>(0,0,0), 0);
	
	
	double g_memsize = sizeof(g[0][0][0])*double(gsizex*gsizex*gsizex) +  sizeof(P*)*p.size();
	double p_memsize = sizeof(p[0])*double(p.size());
	
	all_memsize = (g_memsize+p_memsize)/(1<<20);
	
	fprintf(stderr, "sizeof(P*) = %i [bytes]\n", sizeof(P*));
	fprintf(stderr, "sizeof(vector<P*>) = %i [bytes]\n", sizeof(vector<P*>));
	fprintf(stderr, "sizeof(P) = %i [bytes]\n", sizeof(P));
	
	fprintf(stderr, "pw = %g [Mb]\n", 100.0*sizeof(PW)/(1<<20));
	fprintf(stderr, "g = %g [Mb]\n", g_memsize/(1<<20));
	fprintf(stderr, "p = %g [Mb]\n", p_memsize/(1<<20));
	fprintf(stderr, "mem: %g / %g [Mb]\n", (g_memsize+p_memsize)/(1<<20), (double)memmax);
	
	
	run();


	return 0;
	
}
